package probability;

class ConfluentGemFunction {

	/**
	 * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
	 * 
	 * @since 2.0
	 */
	public static final double GAMMA = 0.577215664901532860606512090082;

	/** Maximum allowed numerical error. */
	private static final double DEFAULT_EPSILON = 10e-15;

	/** Lanczos coefficients */
	private static final double[] LANCZOS = { 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248,
			14.136097974741747174, -0.49191381609762019978, .33994649984811888699e-4, .46523628927048575665e-4,
			-.98374475304879564677e-4, .15808870322491248884e-3, -.21026444172410488319e-3, .21743961811521264320e-3,
			-.16431810653676389022e-3, .84418223983852743293e-4, -.26190838401581408670e-4, .36899182659531622704e-5, };

	/** Avoid repeated computation of log of 2 PI in logGamma */
	private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);

	// limits for switching algorithm in digamma
	/** C limit. */
	private static final double C_LIMIT = 49;

	/** S limit. */
	private static final double S_LIMIT = 1e-5;

	/**
	 * Returns the natural logarithm of the gamma function &#915;(x).
	 * 
	 * The implementation of this method is based on:
	 * <ul>
	 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html"> Gamma Function</a>, equation (28).</li>
	 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html"> Lanczos Approximation</a>, equations (1)
	 * through (5).</li>
	 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on the computation of the convergent
	 * Lanczos complex Gamma approximation </a></li>
	 * </ul>
	 * 
	 * @param x
	 *            the value.
	 * @return log(&#915;(x))
	 */
	public static double logGamma(double x) {
		double ret;

		if (Double.isNaN(x) || (x <= 0.0)) {
			ret = Double.NaN;
		} else {
			double g = 607.0 / 128.0;

			double sum = 0.0;
			for (int i = LANCZOS.length - 1; i > 0; --i) {
				sum = sum + (LANCZOS[i] / (x + i));
			}
			sum = sum + LANCZOS[0];

			double tmp = x + g + .5;
			ret = ((x + .5) * Math.log(tmp)) - tmp + HALF_LOG_2_PI + Math.log(sum / x);
		}

		return ret;
	}

	public static double gamma(double x) {
		return Math.exp(logGamma(x));
	}

	/**
	 * This is from equation 3.5 in the thesis "Computation of Hypergeometric Functions"
	 * 
	 * @param a
	 * @param b
	 * @param z
	 * @return
	 */
	public static double U(double a, double b, double z) {
		return gamma(1 - b) / gamma(1 + a - b) * M(a, b, z) + gamma(b - 1) / gamma(a) * Math.pow(z, 1 - b)
				* M(1 + a - b, 2 - b, z);
	}

	/**
	 * This implementation is from the thesis titled "Computation of Hypergeometric Functions" page 10, method b.
	 * 
	 * @param a
	 * @param b
	 * @param z
	 * @return
	 */
	public static double Mb(double a, double b, double z) {
		double s_1 = 1, s_2 = 1, r = 0, s;
		double TOL = 1E-15;
		s = (a / b) * z;
		for (int j = 2; j < 500; j++) {
			s_1 = s;
			r = (a + j - 1) / (j * (b + j - 1));
			s = s_1 + (s_1 - s_2) * r * z;
			s_2 = s_1;
			if ((Math.abs(s - s_1) / Math.abs(s_1)) < TOL && (Math.abs(s_1 - s_2) / Math.abs(s_2)) < TOL) {
				// System.out.println("break" + Math.abs(s - s_1));
				break;
			}
		}

		return s;
	}

	/**
	 * This implementation is from the thesis titled "Computation of Hypergeometric Functions" page 10, method a.
	 * 
	 * @param a
	 * @param b
	 * @param z
	 * @return
	 */
	public static double M(double a, double b, double z) {
		double aa = 1., s = aa;
		double TOL = 1E-15;
		for (int j = 0; j < 500; j++) {
			aa = aa * (a + (double) j) / (b + (double) j) * (z) / ((double) j + 1.);
			s = s + aa;
			if ((Math.abs(aa) / Math.abs(s)) < TOL) {
				// System.out.println("break" + Math.abs(s - s_1));
				break;
			}
		}

		return s;
	}
}
